Governing equations and semi-Lagrangian framework

In a general point of view, the HySoP library is used to solve continuous systems of the following form:

\[\frac{d}{dt} \displaystyle \int_{\Omega} \mathbf{Q}(x,t) \ d \mathbf{x} = \displaystyle \int_{\Omega} \mathbf{F} (x,t,\mathbf{Q}, \nabla \mathbf{Q}, ...) \ d \mathbf{x}\]

where \(\mathbf{Q}\) denotes the vector of variables and \(\mathbf{F}\) the source term. More precisely, the present library originally lies on the so-called Vortex Methods, which belong to particle (also called Lagrangian) methods. Lagrangian methods differ from Eulerian ones by the fact that the variables \(\mathbf{Q}\) are discretized on a set of particles that follow the dynamic of the system and are displaced with respect to the flow velocity \(\mathbf{u}\). Regarding Vortex Methods, they are used to specifically solve incompressible Navier-Stokes equations in their velocity-vorticity formulation:

\[\begin{split}\begin{align} \frac{\partial\boldsymbol{\omega}}{\partial t} + (\mathbf{u} \cdot \nabla) \boldsymbol{\omega} &= (\boldsymbol{\omega} \cdot \nabla) \mathbf{u} + \dfrac{1}{Re} \Delta \boldsymbol{\omega} + \nabla \times \mathbf{f}_{ext}\\ \Delta \mathbf{u} &= - \nabla \times \boldsymbol{\omega} \end{align}\end{split}\]

The only quantity \(Q\) carried by the particles is the vorticity field \(\boldsymbol{\omega}\), defined in a 3D-Cartesian coordinates system as:

\[\boldsymbol{\omega} = (\omega_x, \omega_y, \omega_z) := \nabla \times \mathbf{u} = \left ( \partial_{y} u_z- \ \partial_{z} u_y\ , \ \partial_{z} u_x - \partial_{x} u_z \ , \ \partial_{x} u_y - \partial_{y} u_x \right )\]

In the above system of governing equations, the first one corresponds to the momentum equation with :

• \((\mathbf{u} \cdot \nabla) \boldsymbol{\omega}\) : the advection term

• \((\boldsymbol{\omega} \cdot \nabla) \mathbf{u}\) : the stretching term (that vanishes in 2D).

• \(\dfrac{1}{Re} \Delta \boldsymbol{\omega}\) : the diffusion term with \(Re\) the Reynolds number.

• \(\nabla \times \mathbf{f}_{ext}\) : the external forcing term that depends on the problem being solved

The second equation, \(\Delta \mathbf{u} = - \nabla \times \boldsymbol{\omega}\), is the Poisson equation allowing to recover the velocity \(\mathbf{u}\) from the vorticity \(\boldsymbol{\omega}\). This equation is derived from the incompressibility condition \(\nabla \cdot \mathbf{u} = 0\) and the definition of the vorticity field \(\boldsymbol{\omega} := \nabla \times \mathbf{u}\).

For a more complete description of the family of models handled by the library, one should rather talk about the resolution of a system of continuous equations consisting of Navier-Stokes equations coupled with \(n\) scalar advection-diffusion equations:

\[\begin{split}\begin{align} \frac{\partial\boldsymbol{\omega}}{\partial t} + (\mathbf{u} \cdot \nabla) \boldsymbol{\omega} &= (\boldsymbol{\omega} \cdot \nabla) \mathbf{u} + \dfrac{1}{Re} \Delta \boldsymbol{\omega} + \nabla \times \mathbf{f}_{ext}\\ \frac{\partial\theta_i}{\partial t} + (\mathbf{u} \cdot \nabla) \theta_i &= \kappa_i \Delta \theta_i \quad \text{for} \ i \in \{1, \cdots, n\}\\ \Delta \mathbf{u} &= - \nabla \times \boldsymbol{\omega} \end{align}\end{split}\]

where \(\kappa_i\) is the constant diffusivity of the scalar \(\theta_i\). In this case, the quantities \(\mathbf{Q}\) carried by the particles are the vorticity field \(\boldsymbol{\omega}\) and the scalar fields \(\theta_i\).

In HySoP, these models are not solved by using a pure Lagrangian approach but rather a semi-Lagrangian method called “remeshed Vortex method” or “remeshed particle method”. Both the momentum equation and the scalar equations can be viewed, at least partially, as advection-diffusion equations, one for the vorticity \(\boldsymbol{\omega}\) and the other for the scalars \(\theta_i\). Those two types of equations can be split into transport and diffusion terms, by relying on so-called operator splitting methods. The idea behind the present numerical method is to split the equations such that each subproblem can be solved by using a dedicated solver based on the most appropriate numerical scheme and by employing a space discretization that is regular enough to be handled by accelerators (GPUs).

Semi-lagrangian (remeshed) particle methods allow to solve advection problems in a Lagrangian way, that is to say directly on particles. In other words the advection of the momentum equation and the scalar advection

\[\frac{\partial\boldsymbol{\omega}}{\partial t} + (\mathbf{u} \cdot \nabla) \boldsymbol{\omega}= 0, \qquad \qquad \qquad \frac{\partial\theta_i}{\partial t} + (\mathbf{u} \cdot \nabla) \theta_i = 0\]

are treated in a Lagrangian way, on each numerical particles \(p\), by solving the following sets of ODEs:

\[\begin{split}\left\{\begin{matrix} \frac{d {\mathbf{x}}_p(t)}{dt} & ={\mathbf{u}}({\mathbf{x}}_p(t), t) \\ \frac{d \boldsymbol{\omega}_p(t)}{dt} &= 0\\ \frac{d \theta^i_p(t)}{dt} &= 0 \end{matrix}\right.\end{split}\]

where the resolution of the first equation updates the numerical particles locations \({\mathbf{x}}_p(t)\) after advection.

Such Lagrangian treatment of the advection equations offers a natural approach, close to the physics, it provides flexible resolution of the non-linear transport problem and ensures stability and low numerical diffusion. It also presents an interesting advantage in terms of computational issues since the Lagrangian advection scheme imposes a CFL stability constraint which is less restrictive than in a Eulerian framework: the Lagrangian CFL condition is indeed based on the velocity gradients and not on a grid size \(\Delta x\), thus allowing the use of larger time steps and also adaptive time steps (\(\Delta t(t)\)).

In order to avoid the distortion of the convected fields, the vorticity and scalar values carried by each particle are distributed (after the advection step) on the neighboring points of an underlying Cartesian mesh. This step is called the “remeshing”. It is done by using remeshing kernels, which are piece-wise polynomial functions, that satisfy desired conservation properties. The vorticity at a node \(i\) of the mesh is thus obtained from the vorticity carried by the neighboring particles \(p\) with weights given by the remeshing kernel \(\Lambda\):

\[\boldsymbol{\omega}_i^{n+1}(x)=\displaystyle \sum_p \boldsymbol{\omega}_p^n(x)\Lambda\left(\frac{x_p^{n+1}-x_i}{\Delta x}\right)\]

In HySoP, the remeshing kernels are denoted \(\Lambda_{m,r}\) where \(r\) corresponds to their regularity \(\mathcal{C}^r\) and \(m\) is the number of preserving moments

Through the projection of the particles on an underlying grid (processed after each advection step) and thank to the operator splitting method, the remeshing process allows the use of eulerian solvers for the treatment of the other operators (ie. stretching, diffusion, external forcing and the Poisson equation). In particular the HySoP library uses Cartesian grids since they are compatible with a wide variety of numerical methods such as finite difference methods (FD) and spectral methods (Fast Fourier Transforms).

In conclusion, the HySoP library is particularly adapted for problems dominated by transport phenomena. However, the operator splitting method on which the library is built allows to handle a wider diversity of problems.